Effects of the Be22W phase formation on hydrogen retention and blistering in mixed Be/W systems
Cao Jin-Li1, He Bing-Ling2, Xiao Wei1, 3, †, Wang Li-Gen3
School of Materials Science and Engineering, University of Science and Technology Beijing, Beijing 100083, China
College of Physics and Electronic Engineering, Xinxiang University, Xinxiang 453003, China
State Key Laboratory of Nonferrous Metals and Processes, General Research Institute for Nonferrous Metals, Beijing 100088, China

 

† Corresponding author. E-mail: wxiao@ustb.edu.cn

Abstract

We have performed first-principles density functional theory calculations to investigate the retention and migration of hydrogen in Be22W, a stable low-W intermetallic compound. The solution energy of interstitial H in Be22W is found to be 0.49 eV lower, while the diffusion barrier, on the other hand, is higher by 0.13 eV compared to those in pure hcp-Be. The higher solubility and lower diffusivity for H atoms make Be22W a potential beneficial secondary phase in hcp-Be to impede the accumulation of H atoms, and hence better resist H blistering. We also find that in Be22W, the attraction between an interstitial H and a beryllium vacancy ranges from 0.34 eV to 1.08 eV, which indicates a weaker trapping for hydrogen than in pure Be. Our calculated results suggest that small size Be22W particles in hcp-Be might serve as the hydrogen trapping centers, hinder hydrogen bubble growth, and improve the resistance to irradiation void swelling, just as dispersed oxide particles in steel do.

1. Introduction

Due to its low atomic number and remarkable neutronic characteristics, beryllium (Be) has been chosen to be the primary plasma facing materials in the international thermonuclear experimental reactor (ITER).[1] As tungsten (W) is always used as the divertor armour materials, the contact of the Be and W layers is unavoidable in nuclear fusion experimental devices. In this case, subsequent re-deposition of a mixed phase Be/W can easily form during physical sputtering and chemical erosion processes, which are induced by deuterium–tritium plasma in the extreme environment of the nuclear reaction.[2] Therefore, H uptake and retention in the mixed material systems are an important issue for reliable extrapolation of in-vessel tritium retention in ITER.

Up to now, various investigations have been devoted to studying the properties of Be/W mixing by experiments,[36] as well as the behavior of hydrogen retention in both pure materials using theoretical simulations.[710] Doerner et al. have found that exposing W to Be-seeded deuterium (D) and D2 plasmas can make the formation of Be–W alloys.[11,12] In addition, the Be2W intermetallic phase can be formed at temperatures above 670 K when Be films are deposited on bulk W and annealed.[13] Conversely, by deposting W films on bulk Be, Wiltner et al. obtained the phase of Be12W within the W diffusion depth at temperatures above 1000 K.[3] Besides the experimental investigations, molecular dynamics (MD) simulations and density-functional theory calculations have also been empolyed to study the hydrogen solubility, diffusion, and interaction with vacancy in pure metals. In the work by Ganchenkova et al., the authors have demonstrated that in bulk Be, hydrogen prefers the basal tetrahedral (BT) interstitial site and the migration energy is 0.40 eV. They have also found that a Be vacancy remains energetically favorable up to five hydrogen atoms.[7] In bulk W, numerous calculations have also indicated that the tetrahedral interstitial site is the favorable one for hydrogen atom and the migration energy is 0.21 eV. In contrast, a W vacancy can trap up to twelve hydrogen atoms.[8] However, the corresponding research on H behavior in Be–W materials is still lacking.[14] From a more fundamental point of view, it is interesting to analyze the influence of the secondary element (W) with rich electrons on the solubility and diffusion of H in the poor-electron element (Be). Furthermore, the physicochemical properties of the mixed materials and alloys are essential input data needed for the modeling of erosion, transport, and deposition in large-scale codes.

Undoubtedly, the studies on H behavior in the mixed materials are supposed to involve many Be/W phases. Based on the investigation of Okammoto et al.,[15] the Be–W phase diagram is very complex, and it includes several stable intermetallic compounds, Be22W, Be12W, and Be2W. Our calculations suggest that the solution energy of W in pure Be is as high as 2.06 eV, and thus if a small amount of tungsten comes into the Be solid, Be22W is supposed to form instead of a solid solution or other intermetallics at the initial stage. With the increase of the tungsten concentration, other stable compounds may finally form. On the other hand, the very low concentration of W, 4.35 at.%, can make the original hcp-structure collapse, which draws our attention to the unexpected interactions between Be and W. In fact, to study H blistering in the mixed Be/W phases, it is necessary to study H behavior in all the stable intermetallic compounds and the corresponding interfaces. However, it will be a big challenge to employ first-principles calculations because of the large mismatch in the constants at the interfaces. As Allouche et al.[14] did for Be12W, here we restrict our study on Be22W, mainly including the effects of Be22W phase formation on hydrogen retention and blistering in the mixed Be/W phase. When high-energy hydrogen ions are implanted into the solid, it inevitably generates Frenkel pairs and cascades of atoms and vacancies, so we will consider H solubility in both interstitial sites and near vacancies. We will organize the remainder of the paper as follows. In Section 2, the theoretical methods and the computational details are described. Section 3 presents and discusses our calculated results of the behavior of the hydrogen retention and diffusion in Be22W. Finally, a short summary is given in Section 4.

2. Methodology

The first-principles density functional theory (DFT) calculations were performed using the Vienna ab initio simulation package (VASP).[1618] The electron–ion interaction was described using the projector augmented wave (PAW) method[19,20] and the exchange and correlation were treated with the generalized gradient approximation (GGA) in the Perdew–Burke–Ernzerhof (PBE) form.[21] The cut-off energy for the plane wave basis set was set to 350 eV for all of the systems. For hcp-Be metals, we employed a 4 × 4 × 3 supercell consisting of 96 atoms, and a 5 × 5 × 5 k-mesh in Monkhorst–Pack scheme for the integration over the Brillouin zone.[22] For cubic Be22W, the calculated supercell contained 184 atoms and the corresponding Monkhorst–Pack mesh was 3 × 3 × 3. The nudged elastic band (NEB) method was adopted to study the hydrogen diffusion behavior.[23] Five images were linearly interpolated between starting and terminal points on the migration track. In the calculations, the supercell parameters were fixed and the optimization of the atomic positions was continued until the total energy of this system was converged to less than 10−4 eV.

3. Results and discussion
3.1. Interstitial H solution and diffusion in bulk Be22W

As shown in Fig. 1, Be22W has a face-centered cubic (fcc) structure (space group No. 227, different from hcp-Be. There are 8 W and 176 Be atoms in one unit cell, in which the W atoms are located at Wyckoff positions 8a, and the Be atoms are at Wyckoff positions 16c, 16d, 48f, and 96g. To clearly display the structure of Be22W, we highlight the W atoms with larger blue spheres, and different Be atoms are represented by blueviolet, fuchsia, violet, and thistle spheres, respectively, from higher symmetry positions to the lower ones. Compared to the unit cell of Be, we can find that with the addition of 4.35 at.% W atoms, the original hcp structure collapses, and the Be atoms are distributed into four types of lattices with different bonding environments. The lattice parameter of fcc Be22W is 11.56 Å, and its formation energy with reference to pure metals is −0.04 eV/atom, which testifies that the compound can be easily formed. These results are in good agreement with the previous DFT calculations.[24]

Fig. 1. (color online) The unit cell of Be22W: W8a, Be16c, Be16d, Be48f, and Be96g atoms are in blue, blueviolet, fuchsia, violet, and thistle spheres, respectively. Various interstitial (Int) sites are denoted by green spheres with letters A–G.

At first, we focus on the hydrogen solution at the interstitial sites of the Be22W matrix. Seven different interstitial sites have been considered in the calculations. These sites are denoted by green spheres with the letters A–G in Fig. 1 and also zoomed in Figs. 2(a)2(g) for clarity. We define the solution heat of interstitial hydrogen as the energy needed to embed it into the interstitial position from the free gas state

where E(Be176W8H) and E(Be176W8) are the total free energies of systems with and without an interstitial H atom, respectively, is the chemical potential for H, and equals to one half of the free energy of the isolated H2 molecular, −3.38 eV, which is −3.25 eV with zero-point energy correction. At very low electron density zone, the energy cost to embedding an H atom into an electron gas is in proportion to the electron density, and H always prefers the low electron density space.[25] However, as an interstitial atom in a solid, the factor of volume and the interaction with surrounding metal atoms are also very significant. From Figs. 2(a)2(c), we can see that the interstitial sites A and B are regular octahedral, while the site C is irregular. In addition, the site A is embraced by six chemically identical Be atoms, and the Be atoms around the sites B and C are distinct. The polyhedral volumes are 5.07 Å3, 5.31 Å3, and 5.83 Å3 for the interstitial sites A, B, and C, respectively. For other interstitials, the discrepancy of the bonding environment is more obvious. The tetrahedral interstice D in Fig. 2(d) is surrounded by one W atom and three Be atoms with a volume of 1.52 Å3. The interstitial site E in Fig. 2(e) has 2 nearest neighbors (NNs), Be and W atoms with a distance of 1.25 Å, and 6 next-nearest-neighboring (NNN) Be atoms with a distance of 2.16 Å. The interstitial sites F in Fig. 2(f) and G in Fig. 2(g) in a Be–W or Be–Be triangle with very short bond lengths are also included in our work due to the fact that the BT interstitial site possesses the lowest solution heat of 1.47 eV for hydrogen in hcp-Be.

Fig. 2. (color online) (a)–(g) Zoom in to see different interstitial sites A–G in Fig. 1. Note that only the local structures, bonds, and bond lengths around the defect are displayed. For an explanation of the color spheres representing different atoms, see Fig. 1.

According to formula (1), we calculate the hydrogen solution energies at the interstitial sites, and list the results in Table 1. For the interstitial site E, it has a high solution energy of 1.74 eV because of severe compression induced by the short distance between H and its NN metal atoms. The solution energies at the interstitial sites A and G are almost equal, namely, 1.54 eV and 1.56 eV. In other interstices, it ranges from 0.99 eV to 1.39 eV, lower than that in hcp-Be, which has been confirmed by the previous calculations.[26] Owing to the low atomic mass, zero-point energy (ZPE) correction of H should also be taken into account. In Table 1, we list the corresponding solution energies at different interstices corrected with ZPE and find that the trend cannot be changed. Thus we can confirm that hydrogen prefers to occupy the interstitial site C followed by D, B, and F in Be22W, in which the lowest solution energy is 0.49 eV smaller than that in hcp-Be.

Table 1.

The solution energies ( and corresponding values ( considering zero point energies of interstitial hydrogen in Be22W. The different interstitial sites are schematized by green spheres in Figs. 1 and 2.

.

As we all know, the concentration of point defects such as interstitial H in the thermal equilibrium is determined by the formula

where A is the coordination number, is the solution energy of the interstitial H,[9,26] and is the Boltzmann constant. Thus, the study of the early-stage growth of H clusters can be performed by thermodynamic analysis. In Fig. 3, we directly compare the H concentrations in Be22W and hcp Be. Due to different lattice structures of Be22W and Be, here we consider their concentration in unit volume. Since there are six interstitial BT per Be unit cell with a volume of 47.52 Å3, and 96 interstitial C in the calculated Be22W cell with a volume of 1544.80 Å3, A in formula (2) is 0.13 Å−3 for hcp-Be, and 0.06 Å−3 for Be22W. We depict the H concentration in the denary logarithm as a function of temperature in Fig. 3. It shows that the concentration of interstitial H increases with the elevated temperature, and in Be22W it is remarkably higher than that in pure Be, for instance, the concentration of H is 8.33 × 10−11−3 in Be22W, 1.26 × 10−14−3 in Be at 600 K, and 7.66 × 10−8−3 in Be22W, 2.71 × 10−10−3 in Be at 900 K, which are labeled by rectangles in Fig. 3. It means that the occupation probabilities of H in Be22W can reach up to hundreds or thousands of times higher than those in a Be matrix. A high local concentration of interstitial H can appear at the formed Be22W compound in the Be matrix. Therefore, these inherent interstitial sites in Be22W can act as traps for hydrogen.

Fig. 3. (color online) Interstitial hydrogen concentration in denary logarithm as a function of temperature in different solids, Be22W (blue) and Be (red). The points at 600 K and 900 K in both solids are labeled by rectangles.

In order to detailedly elucidate the effect of tungsten addition on the H behavior in Be22W, we should also consider the kinetics of H accumulation inside the Be22W solid, i.e., the diffusion and migration paths of hydrogen. Here we only consider the migration paths and barrier between the relatively stable interstitial sites using the NEB method. According to the solution energies obtained above, the initial and final positions for the diffusion paths are both interstitial sites C in Be22W.

As shown in Fig. 4(a), the interstitial sites C in the unit cell are depicted in green spheres, which form a regular hexagon plotted by green lines in the vicinity of W atoms. Each W atom has four interstice C-formed hexagons, and the center of each hexagon is Be16c. Therefore, H can migrate along the regular hexagon between two NN interstitial sites C with a distance of 1.73 Å. The migration energy profiles for this path () are described by blue circles on the left in Fig. 4(b). It shows that the H atom has to overcome an energy barrier of 0.50 eV when it migrates from C1 to C2 positions.

Fig. 4. (color online) Migration of interstitial hydrogen: (a) the NN interstitial sites C (green spheres) form a regular hexagon, located in the vicinity of W atoms (blue spheres). Interstitial sites G (yellow spheres) are located in the middle of the NNN interstitial site C. Note that we only represent W atoms by blue spheres for observation convenience. (b) Migration energy profiles for hydrogen migration in Be22W are estimated using NEB calculation. The C2 site is the NN site of interstitial site C1 with a distance of 1.73 Å, and C3 is the NNN site of site C2 with a distance of 2.04 Å.

Furthermore, H needs to diffuse in full migration steps not only along the regular hexagon. By full we mean the center of mass of an H atom moves from position (an interstitial site C2, for instance) to , where , , and are the unit vectors with a length of lattice constant. Then, we should consider another kind of hydrogen diffusion and migration path between two NNN interstitial C sites, i.e., from C2 to C3. In Fig. 4, we show that the migration from C2 to C3 is 2.04 Å far away, and during the migration H shall pass through an interstitial site G depicted by the yellow sphere, which is just located in the center of two sites. As displayed in Fig. 4(b), the migration energy barrier along this path is 0.52 eV, which is indeed determined by the solution energy difference between C and G positions. Through the calculations of two types of diffusion and migration paths, we can obtain a migration energy barrier of 0.52 eV for H inside Be22W.

As a comparison, we calculate the migration energy barrier of H inside the Be matrix, which is 0.39 eV from BT to nearby BT through octahedral (O) sites. The result is consistent with the previous theoretical values (0.38 eV,[7] 0.40 eV,[26] 0.41 eV[27]). This means that the migration energy barrier inside Be22W is 0.13 eV higher than that in the Be matrix. This lower diffusivity for H atoms in Be22W can impede the movement of H atoms.

3.2. H trapping of vacancy in Be22W

When high-energy hydrogen ions are implanted into the solid, it will generate Frenkel pairs, cascades of atoms, and vacancies in the materials. Vacancies as a very low electron density zone always act as traps for hydrogen. Therefore, the studies of the vacancy formation and its interaction with H in Be22W are significant.

At the very beginning, we consider different types of monovacancy by removing one W atom or one of four chemically distinct Be atoms respectively in Be22W, and their formation energies can be defined as

or
where E(Be176W7) and E(Be175W8) are the total free energies of the cell with one W and Be vacancy, and and are the chemical potentials for W and Be, respectively. As Be22W in our work is formed at the interface between pure beryllium and tungsten solids, the chemical potential of W and Be has been taken as its value in bcc-W and hcp-Be solid.

Using formulas (3) and (4), we calculate the formation energy of monovacancy for tungsten and beryllium in Be22W, and list these results in Table 2. In contrast with pure crystals, the formation energies of a vacancy in Be and W are also included. Since W has a much larger cohesive energy than Be, the formation of a vacancy in W can be more difficult. The calculated values are 3.26 eV and 1.09 eV, respectively, for W and Be, consistent with the previous DFT calculations[9,10,14] and the experimental values.[28] For the W monovacancy in Be22W, it turns out that the formation energy is 4.91 eV, which shows a remarkable increase compared to the case of pure W. For the Be vacancy in Be22W, the results are a bit more complicated. Four different types of vacancies can be formed by removing one Be16c, Be16d, Be48f, and Be96g atom respectively in the supercell. Our calculated results show that the formation energies of Be16c and Be16d vacancies increase to 1.88 eV and 1.26 eV, while the Be48f and Be96g ones decrease to 0.86 eV and 1.03 eV in comparison with the pure Be. Generally, Be16c and Be16d vacancies in Be22W can be relatively hard produced, which agrees with the work in Be12W by Allouche et al. However, the change trends of Be48f and Be96g vacancies are different because of the different bonding environments. In Table 2, it shows in Be12W, three chemically distinct Be vacancies all have higher formation energies than those in pure Be, while the W vacancies have lower formation energies.[14]

Table 2.

Formation energies of monovacancy in Be22W. For comparison, the values in Be12W, pure W, and pure Be are also listed.

.

To gain deeper insight into the change trend of vacancy formation in Be22W, we now turn to the electronic structure analyses. We plot the electronic densities of states (DOS) of W atoms in Be22W and pure W in Fig. 5(a), and Be atoms in Be22W and pure Be metal in Fig. 5(b). Compared with the pure W, the broadened bands including conduction and valence bands in the DOS of tungsten in Be22W indicate that the W atoms suffer the compression. This is consistent with the atomic structure, in which the NN distances are 2.50 Å in Be22W, shorter than those in pure W, 2.75 Å. Thus due to the release of the elastic energy, the tungsten vacancy should be easier to form. However, we can also see an additional peak at −0.6 eV marked by the dot ellipse in Fig. 5(a), which suggests that there are interactions between W and Be atoms in Be22W. Furthermore, the DOS at the Fermi energy level has a decrease, see the inset in Fig. 5(a). As a consequence, the formation of tungsten vacancy in Be22W becomes harder. In Fig. 5(b), the p-state peaks of the valence band for Be16c and Be16d in the range [−4.00, −1.60] eV, which are also d-state peaks of the valence band of W in Be22W, have a remarkable increase. Especially for Be16c, the increase of peaks at −2.4 eV is simultaneous with W, and this declares the stronger interaction between them. This is mainly determined by the atomic structure, tungsten atoms being the NN of Be16c with a short distance of 2.50 Å. For Be48f and Be96g, the top of the valence band has a weak increase, and resembles pure Be, thus the formation energy of vacancies is close to the pure metal. In addition, the DOSs at the Fermi energy level of Be16c and Be16d are lower than those of Be48f and Be96g, which suggests that Be16c and Be16d are harder to form.

Fig. 5. (color online) Comparison of electronic density of states of (a) W atoms in Be22W and pure W, (b) Be atoms in Be22W and pure Be. The Fermi energy is set to zero denoted by black dash lines, and the DOS nearby the Fermi energy level ranging from −1 eV to 1 eV are drawn in the inset.

Hydrogen can facilitate the preservation of superabundant vacancies in a large variety of metals and alloys,[29] and vice versa, vacancies have been found to trap hydrogen.[8,9] Therefore, it is necessary to investigate the interaction between H and the vacancy. We define the trapping energy ( to characterize the energy required for moving an H atom into the vicinity of a beryllium vacancy from a remote interstitial site C

where is the total free energy of systems when H is near one Be atom vacancy. Here, we just consider a beryllium monovacancy due to the large formation energy of the tungsten vacancy. In the alloy calculations, the H atom is placed at and off the center of the monovacancy, and the optimized atomic configurations suggest that H atoms prefer to stay at 1.71 Å, 0.98 Å, 1.17 Å, and 0.92 Å away from the center of Be16c, Be16d, Be48f, and Be96g vacancy, respectively. We list the distances between H and the vacancy and the calculated trapping energies of various beryllium vacancies in Table 3. The trapping energy of Be96g vacancy is 1.08 eV, slightly lower than that of pure beryllium, 1.14 eV. Nevertheless, for three other types of beryllium vacancies, the trapping energies of beryllium vacancy, Be16c, Be16d, and Be48f for H atoms are 0.34 eV, 0.39 eV, and 0.63 eV, respectively, markedly lower than that of pure beryllium. Among them, the low trapping energies of Be16C and Be16d monovacancy can be attributed to their large vacancy formation energies. The weaker trapping of a vacancy for H in Be22W provides a less driving force to form H clusters than that in pure Be solids.

Table 3.

Calculated distances between hydrogen and vacancy ( and trapping energies of various beryllium vacancies in Be22W.

.
4. Conclusion

Using first-principles density calculations, we have systematically investigated the distribution and migration of H inside Be22W solid. The solution energy of interstitial H atoms in Be22W is 1.06 eV, which is 0.49 eV smaller than that in hcp-Be. The migration barrier of interstitial H is 0.52 eV from the most stable interstitial site C to its NNN site, larger than that in Be solid. Therefore, compared to H in pure Be solid, the larger solubility and slower diffusion of H atoms in Be22W solid make it a potential beneficial secondary phase in hcp-Be to impede hydrogen diffusion. On the other hand, the interaction of beryllium and tungsten in Be22W increases the formation energy of W, Be16c and Be16d vacancies and decreases it for Be48f and Be96g in comparison with the pure Be. Furthermore, we have obtained the hydrogen trapping energies of distinct beryllium vacancies in Be22W, including Be16c, Be16d, Be48f, and Be96g, which are 0.34 eV, 0.39 eV, 0.63 eV, and 1.08 eV. All are lower than that in pure Be, 1.14 eV, which provides a less driving force for the accumulation and clustering of H atoms than that in pure Be solids.

Our calculations suggest that interstitial H atoms have a lower solution energy, a higher migration energy, and are harder to accumulate into the vacancy in Be22W than the matrix. Based on our calculations, H prefers Be22W over pure Be solid, which can slow the growth of H clusters and bubbles in Be solid. Thus, H bubbles are dispersed and the concentrated stress may be released below the limited strength of the materials. Consequently, H embrittlement or hydrogen-related irradiation swelling could be suppressed. We can hope that the presence of Be22W is potentially very effective to serve as a dispersed trapping center to retard the growth of H bubbles, just as the suppressing effect of oxide particles on the growth of He and H bubble in ODS steels,[30,31] and enhancing the resistance to irradiation void swelling and H embrittlement.

Reference
[1] Merola M 2010 Fusion Eng. Des. 85 2312
[2] Perkins F W 1999 Nucl. Fusion 39 2137
[3] Wiltner A Kost F Lindig S Linsmeier Ch 2007 Phys. Scr. T128 133
[4] Linsmeier Ch 2007 J. Nucl. Mater. 363�?65 1129
[5] Baldwin M J 2007 J. Nucl. Mater. 363�?65 1179
[6] Lasa A Heinola K Nordlund K 2014 Nucl. Fusion 54 083001
[7] Ganchenkova M G Borodin V A Nieminen R M 2009 Phys. Rev. 79 134101
[8] Ohsawa K Goto J Yamakami M Yagi M 2010 Phys. Rev. 82 184117
[9] Zhang P B Zhao J J Wen B 2012 J. Phy.: Condens. Matter. 24 095004
[10] Xiao W Luo G N Geng W T 2012 J. Nucl. Mater. 421 176
[11] Doerner R P Baldwin M J Causey R A 2005 J. Nucl. Mater. 342 63
[12] Doerner R P 2007 Phys. Scr. T128 115
[13] Wiltner A Linsmeier C 2005 J. Nucl. Mater. 337 951
[14] Allouche A Fernandez N Ferro Y 2014 J. Phys.: Conden. Matter. 26 315012
[15] Okammoto H Tanner L E 1991 Phase Diagrams of Binary Tungsten Alloys Calcutta Indian Institute of Metals
[16] Kresse G Furthmüller J 1996 Phys. Rev. 54 11169
[17] Kresse G Hafner J 1993 Phys. Rev. 47 558
[18] Kresse G Furthmüller J 1996 Comp. Mater. Sci. 6 15
[19] Kresse G Joubert D 1999 Phys. Rev. 59 1758
[20] Blöchl P E 1994 Phys. Rev. 50 17953
[21] Perdew J P Burke K Ernzerhof M 1996 Phys. Rev. Lett. 77 3865
[22] Monkhorst H J Pack J D 1976 Phys. Rev. 13 5188
[23] Henkelman G Uberuaga B P Jónsson H 2000 J. Chem. Phys. 113 9901
[24] Jain A 2013 APL Mater. 1 011002
[25] Puska M J Nieminen R M Manninen M 1981 Phys. Rev. 24 3037
[26] Zhang P B Zhao J J Wen B 2012 J. Nucl. Mater. 423 164
[27] Allouche A Oberkofler M Reinelt M Linsmeier Ch 2010 J. Phys. Chem. 114 3588
[28] Baskes M I Johnson R A 1994 Model. Simul. Mater. Sci. 2 147
[29] Fukai Y 2003 Phys. Scr. T103 11
[30] Hsiung L L 2010 Phys. Rev. 82 184103
[31] Esteban G A 2007 Fusion Eng. Des. 82 2634